$natural_history
to
from Healthy CVD Dead
Healthy 0.99956 0.00041087 0.000027397
CVD 0.00000 0.99972606 0.000273935
Dead 0.00000 0.00000000 1.000000000
In the presence of competing events (mortality), we obtain different results for different cycle lengths! Notice how the state occupancy in CVD differs in the right panel of the figure below.
$natural_history
to
from Healthy CVD CVDDeath Dead
Healthy 0.12797 0.77687 0.00000 0.095163
CVD 0.00000 0.36788 0.63212 0.000000
CVDDeath 0.00000 0.00000 1.00000 0.000000
Dead 0.00000 0.00000 0.00000 1.000000
Discretizing a Continuous Process
To accurately capture the underlying continuous time process, there should be a transition probability defined for a seemingly impossible transition: Healthy CVD Death.
$natural_history
to
from Healthy CVD CVDDeath Dead
Healthy 0.12797 0.77687 0.00000 0.095163
CVD 0.00000 0.36788 0.63212 0.000000
CVDDeath 0.00000 0.00000 1.00000 0.000000
Dead 0.00000 0.00000 0.00000 1.000000
Discretizing a Continuous Process
Discretizing a Continuous Process
Why Don’t the Conversion Formulas Work?
When we use the standard rate-to-probability conversion formula (i.e., \(p_{H,CVD}= 1 - \exp(-r_{H,CVD} t)\) ) we obtain the marginal transition probability.
This marginal probability reflects the union of all possible transitions that start Healthy and then transition to and from the CVD state in the cycle.
Why Don’t the Conversion Formulas Work?
That is, \(p_{H,CVD}= 1 - \exp(-r_{H,CVD} t)\) captures the probability of following scenarios occurring in a cycle:
Individual transitions from healthy to CVD.
Individual transitions from healthy to CVD to CVD Death.
Individual dies from background causes before they would have developed CVD in the same cycle.
$natural_history
to
from Healthy CVD CVDDeath Dead
Healthy 0.12797 0.77687 0.00000 0.095163
CVD 0.00000 0.36788 0.63212 0.000000
CVDDeath 0.00000 0.00000 1.00000 0.000000
Dead 0.00000 0.00000 0.00000 1.000000
But we actually need to distribute out the marginal probability so that we’re left only with the probability that someone transitions to CVD and remains there at the end of the cycle.
As constructed, we have ruled out the possibility of compound transitions (Healthy CVD CVDDeath) in a single cycle.
$natural_history
to
from Healthy CVD CVDDeath Dead
Healthy 0.12797 0.77687 0.00000 0.095163
CVD 0.00000 0.36788 0.63212 0.000000
CVDDeath 0.00000 0.00000 1.00000 0.000000
Dead 0.00000 0.00000 0.00000 1.000000
This results in overcounting CVD state occupancy in some cases (people who develop CVD and die soon after) and undercounting CVD state occupancy in others (people who die of background causes before they would have developed CVD in the same cycle).
Discretizing a Continuous Process
Changing to shorter cycle does not get us out of this problem.
Quantified Absolute Error
Error after one cycle
truth <-Vectorize(function(r1, r2, t){# Embed properly x <-matrix(c(-r1-r2, r1, r2, rep(0, 6)), nrow=3, ncol=3, byrow=TRUE) (c(1, 0, 0) %*%expm(x*t))[,2:3]})misspecify <-Vectorize(function(r1, r2, t){# Common Advice# RL Fleurence, Christopher S. Hollenbeak. Rates and Probabilities in Economics, Pharmacoeconomics 2007; 25 (1) 3-6 tp1 <-1-exp(-r1*t) tp2 <-1-exp(-r2*t)# THIS CONVERSION CAN RESULT IN IMPOSSIBLE VALUES!if(tp1 + tp2 >1) return(NA)# Mispecify x <-matrix(c(1-tp1-tp2, tp1, tp2, 0, 1, 0, 0, 0, 1), nrow=3, ncol=3, byrow=TRUE) (c(1, 0, 0) %*% x)[,2:3]})par(mfrow=c(1,1))abs.error <-Vectorize(function(r1, r2, t=1){misspecify(r1, r2, t) -truth(r1, r2, t)})x <-seq(0.01, 0.7, by =0.01)y <- xz <-outer(x, y, FUN=function(x, y) {abs.error(x, y)[1,]})contour(x,y,z, xlab="Rate 1 A -> B", ylab="Rate 2 A -> C",main="Misspecified Competing Rates",sub="Absolute Error Contour State B, t=1")
Quantified Relative Error
rel.error <-Vectorize(function(r1, r2, t=1){ x <-truth(r1, r2, t)100*(misspecify(r1, r2, t) - x)/x})par(mfrow=c(1,2))z <-outer(x, y, FUN=function(x, y) {rel.error(x, y)[1,]})contour(x,y,z, xlab="Rate 1 A -> B", ylab="Rate 2 A -> C",main="",sub="Relative Error Contour (%) State B, t=1")z <-outer(x, y, FUN=function(x, y) {rel.error(x, y)[2,]})contour(x,y,z, xlab="Rate 1 A -> B", ylab="Rate 2 A -> C",main="",sub="Relative Error Contour (%) State C, t=1")par(mfrow=c(1,1))title(main="Misspecified Competing Rates")
For rates < 0.1 of the chosen time step the error is <5%.
This error will accumulate over each timestep.
Comparisons Across Strategies Bail You Out
Common CEA approaches compare strategies, i.e., \(f(\theta_A) - f(\theta_{B})\)
Net monetary benefit.
Net health benefit.
Incremental cost effectiveness ratio.
With error it becomes:
Comparisons Across Strategies Bail You Out
Common CEA approaches compare strategies, i.e., \(f(\theta_A) - f(\theta_{B})\)
If \(\epsilon_A - \epsilon_{B} \sim 0\) the impact to the decision threshold is minimal.
These decision thresholds are semi-robust against this misspecification.
Restated: Single model runs are biased under misspecification, but comparison across strategies helps cancel this out.
Rate Misspecification Conclusions
In the presence of competing events (mortality), using standard conversion formulas will yield incorrect transition probability matrices.
Why Does This Matter?
Without properly embedding, you are no longer modeling a continuous time process.
So your model results will differ compared to a continuous time model (e.g., discrete event simulation)
Microsimulation suffers the same problem because the patient-level transition probability matrix is embedded in a time step and requires Markov Rate conversion.
Rate Misspecification Conclusions
Existing publication’s decision thresholds are generally okay if they misspecified rates.
Check the maximum rate of an event as a rough estimate of potential bias.
Rule of thumb, on the error plot use the estimated rate for full time scale of simulation, i.e. one step is entire simulation.
Rate Misspecification Error Decision Threshold
The amount of error for a decision threshold is the difference in contour. In this example: \(0.07 - 0.02 = 0.05\).
truth <-Vectorize(function(r1, r2, t){# Embed properly x <-matrix(c(-r1-r2, r1, r2, rep(0, 6)), nrow=3, ncol=3, byrow=TRUE) (c(1, 0, 0) %*%expm(x*t))[,2:3]})misspecify <-Vectorize(function(r1, r2, t){# Common Advice# RL Fleurence, Christopher S. Hollenbeak. Rates and Probabilities in Economics, Pharmacoeconomics 2007; 25 (1) 3-6 tp1 <-1-exp(-r1*t) tp2 <-1-exp(-r2*t)# THIS CONVERSION CAN RESULT IN IMPOSSIBLE VALUES!if(tp1 + tp2 >1) return(NA)# Mispecify x <-matrix(c(1-tp1-tp2, tp1, tp2, 0, 1, 0, 0, 0, 1), nrow=3, ncol=3, byrow=TRUE) (c(1, 0, 0) %*% x)[,2:3]})par(mfrow=c(1,1))abs.error <-Vectorize(function(r1, r2, t=1){misspecify(r1, r2, t) -truth(r1, r2, t)})x <-seq(0.01, 0.7, by =0.01)y <- xz <-outer(x, y, FUN=function(x, y) {abs.error(x, y)[1,]})contour(x,y,z, xlab="Rate 1 A -> B", ylab="Rate 2 A -> C",main="Misspecified Competing Rates",sub="Absolute Error Contour State B, t=1")points(x[8], y[63], col='red', lwd=4, cex=2)points(x[63], y[37], col='red', lwd=4, cex=2)lines(c(x[63], x[63]), c(y[37], y[9]), lty=2, lwd=2)